{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pyccl as ccl\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up cosmologies in CCL\n",
    "p1 = ccl.Parameters(Omega_c=0.25, Omega_b=0.05, h=0.7, sigma8=0.80, n_s=0.96, w0=-1.0, wa=0.0)\n",
    "p2 = ccl.Parameters(Omega_c=0.25, Omega_b=0.05, h=0.7, sigma8=0.80, n_s=0.96, w0=-0.9, wa=0.0)\n",
    "p3 = ccl.Parameters(Omega_c=0.25, Omega_b=0.05, h=0.7, sigma8=0.80, n_s=0.96, w0=-0.9, wa=0.1)\n",
    "p4 = ccl.Parameters(Omega_c=0.25, Omega_b=0.05, h=0.7, sigma8=0.80, n_s=0.96, w0=-0.9, wa=0.1)\n",
    "p5 = ccl.Parameters(Omega_c=0.25, Omega_b=0.05, h=0.7, sigma8=0.80, n_s=0.96, w0=-0.9, wa=0.1)\n",
    "\n",
    "p1.parameters.Omega_g = 0\n",
    "p2.parameters.Omega_g = 0\n",
    "p3.parameters.Omega_g = 0\n",
    "p4.parameters.Omega_g = 0\n",
    "p5.parameters.Omega_g = 0\n",
    "\n",
    "p4.parameters.Omega_l = 0.65\n",
    "p5.parameters.Omega_l = 0.75\n",
    "\n",
    "# define the cosmology model to use\n",
    "p = p3\n",
    "cosmo = ccl.Cosmology(p)\n",
    "a = 1\n",
    "beta = 0.5102500557881973\n",
    "rsd_bench = np.loadtxt('../tests/benchmark/model3_xiRSD.txt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEeCAYAAACg8JNZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XmYXGWd9vHvL/vWSQjdhE0SQkgT9iVsChpQFJQgg+wOgjKACs68wKggjO/MOwiiIq7IACNRERAEImAUXAiLCUhQWdIRSAKSELKQfSXpzu/946lDVVdXdVfVqe1U35/rqquqT52qevp0dd31rMfcHRERkVL1qXUBREQk2RQkIiISi4JERERiUZCIiEgsChIREYlFQSIiIrEoSKRhmdlkM1sU4/E3m9l/lLNMeV5nqpldk7p9tJm9nHFfq5n91czWmdm/mtlgM3vIzNaY2b2VLptIIfrVugAi3TGz14HRQAewHvgtcIm7ry/z65wH/Iu7HxVtc/fPlvM1CuHuTwKtGZu+BMxw94MAzOwcwvHY3t3bq10+kVxUI5EkmOLuw4ADgYOAK2tcnmoaA8zJ+vmVUkLEzPTFUSpCQSKJ4e5LgEcIgQKAmQ00s2+Z2RtmtjTVHDU41+PN7Aozm59qJmozs39KbZ8I3AwcaWbrzWx1antmk9NcMzsx47n6mdnbZnZw6ucjzGymma02s+fNbHK+38PMDjKzv6TK8QtgUMZ97zbHmdkfgWOAH6TKdRfwVeCM1M/np/b7TKp8q8zsETMbk/F8bmYXm9mrwKupbXuZ2e/MbKWZvWxmp2fsP9XMfmhmv06V7xkz2yPj/n0yHrvUzL6S2t4n4/iuMLN7zGxUN39OaSAKEkkMM9sVOAGYl7H5emACIVzGA7sQPmxzmQ8cDYwA/gu4w8x2cve5wGeBWe4+zN1H5njsXcBZGT9/BHjb3f9iZrsAvwauAUYB/w7cZ2YtOX6HAcA04Gepfe8FPpGrsO5+LPAkoSlvmLufBVwL/CL18/+a2cnAV4BTgJbU/ndlPdXJwOHA3mY2FPgdcCewQ+p3usnM9snY/6zU8dmOcKy/lip7E/B7QvPizoTj/YfUY/419TofSN23Cvhhrt9LGo+CRJJgmpmtAxYCy4D/C2BmBlwAXOruK919HeGD9sxcT+Lu97r7Ynff5u6/IHxDP6zAMtwJnGRmQ1I/n53aBvDPwHR3n5567t8Bs4GP5nieI4D+wHfcfau7/xJ4tsAy5HIRcJ27z001d10LHJhZK0ndv9LdNwEnAq+7++3u3u7ufwHuA07N2P9+d/9z6vl+TroGeCKwxN1vcPfN7r7O3Z/JKMdV7r7I3d8B/hM4Vc1pvYOCRJLgZHdvAiYDewHNqe0twBDguVST0mrCt+UuNQEAM/uUmf0tY999M56rW+4+D5gLTEmFyUmkg2QMcFr0vKnnPgrYKcdT7Qy86Z1XS/1HIWXIYwzw3YzXXQkYoWYWWZi1/+FZZf0ksGPGPksybm8EhqVuv4dQq8tXjgcynnMuYYDE6BJ/L0kQfVuQxHD3x81sKvAtQjPK28AmYB93f7O7x6a+od8KfJDQhNVhZn8jfOgCFLIMdtS81QdoS4ULhA/qn7n7BQU8x1vALmZmGWGyG/k/oHuyEPiau/+8m30yf7eFwOPuflyJr3VWN/d9xt3/VMLzSsKpRiJJ8x3gODM70N23EcLhRjPbAcDMdjGzj+R43FDCB+ry1H6fJtRIIkuBXVN9GPncDXwY+Bzp2gjAHYSaykfMrK+ZDUp1mu+a4zlmAe3Av6Y67E+h8Oa1XG4Groz6OMxshJmd1s3+DwMTzOwcM+ufuhyaGnDQk4eBHc3s/6QGOTSZ2eEZ5fha1KRmZi1m9vEYv5ckiIJEEsXdlwM/BaKJgl8mdAg/bWZrCZ3BrTke1wbcQPggXwrsB2R+e/4jYZjtEjN7O89rv5V6/HuBX2RsXwh8nNDpvZzw7fyL5Pj/cvcthI7x8wgd0mcA9xfyu+cp0wOEAQd3p37/lwgDEvLtv44QhmcCiwnNWNcDAwt4rXXAccCU1ONeJYwqA/gu8CDwaKo/62lCB7/0AqYTW4mISByqkYiISCwKEhERiUVBIiIisShIREQkFgWJiIjE0tATEs1sCjClqanpggkTJpT0HBs2bGDo0KHlLZh0oeNceTrG1dFIx/m55557291zrhSRqVcM/500aZLPnj27pMfOmDGDyZMnl7dA0oWOc+XpGFdHIx1nM3vO3Sf1tJ+atkREJBYFiYiIxKIgERGRWBo6SMxsipndsmbNmloXRUSkYTV0kLj7Q+5+4YgRI2pdFBGRhtXQQSIiIpWnIBERqXfr18P8Us99VnkKEhGRevfBD8L48bUuRV4KEhGRevfnP9e6BN1q6CDRqC0RaSh1uhJJQweJRm2JSEPZtq3WJcipoYNERKShtLfXugQ5KUhERJJi69ZalyAnBYmISFKoRiIiIrEoSEREJBY1bYmISCyqkVSf5pGISENRkFSf5pGISENR05aIiMSiGomIiBQtc1kUBYmIiBRt8+b0bTVtiYhI0TZtSt9WjURERIqmIBERkVgyg0RNWyIiUrTMPhLVSKpPExJFJPEyayEKkurThEQRSbzM8FDTloiIFC0zSFQjERGRoilIREQkFjVtiYhILOpsFxGRWNS0JSIisahpS0REYklAjaRfrQtQz6ZOhQcfHM/MmTB6NOy4Y/p6hx1gwIBal1BEGp6CJNn+9jd49NEdeeCB3PdfdBHcfHN1yyQivUwCmrYUJN34znfg5JOf4ogjJrN0KSxZEi5Ll8IPfwhPPVXrEopIw0vAqC0FSQEGDYIxY8Il8txzMG1a7cokIr1EZni8/HLtytENdbaXqKUFVqyAbdtqXRIRaWiZQTJ1Kjz2WM2Kkk9DB0klV/9taYGODli1quxPLSKSlt2ctWBBYY/70pfgjjvKX54cGjpIKrn6b0tLuF6+vOxPLSKSlh0kw4b1/JjXX4dvfhPOOQdeeaUixcrU0EFSSQoSEamK7JFaGzf2/Jj77gvXe+wRmk4qTJ3tJVKQiEhVZNdI1q3r+TFz5sDOO8O8eZUpUxbVSEqkIBGRqigmSDZuBPfQtDV2bCVL1YmCpETNzeFaQSIiFRUFyfTp4TpXkKxdC1dcAUOHwsiRYWRX5nyFClOQlGjgQBg+XEEiIhXW3g59+8IJJ8CIEXD99XDNNZ33+Y//CNshhArAe95TtSIqSGJoaSk+SNaurdvJqSJSj9rboV+qOzuqjVx3Xfr+V1+Fu+7q+jj3ypctRUESQylBst9+cMMNlSmPiDSgrVvTQRLNgB40KH3/xIm5P4guu6zyZUtRkMRQbJCsXw9vvFH4fCIREdrboX//ztuiIFm2LPfw3jvuCMuUV4mCJIZig2TJknBdyOg9ERGgc9NWJAqS0aNzP2bnnStbpiwKkhiiICm0KfKtt8K1gkRECpYrSBYsgA98IP9jtt++smXKoiCJoaUlNF9GgyR6EtVICt1fRKRTkGR2qj/xRO79J06EPfesfLkyKEhiKHZSomokIlK0zM72M8+ET30q/74PPghtbTB4cHXKlqIgiaHYIFEfiYgULbtp6+KLYa+94N//vfN+11wDU6ZUt2wpCpIYSg0SNW2JSMGyR20ddhjMnQvf+Eb4MDnvvDB666qralZELdoYg5q2RKTicnW2A5hBUxPcfnv1y5RFNZIYSq2RbNpUvtnt//gHvPe9oZY7e3ZVJ7OKSDXkC5I6oiCJYciQcCm2RgLlq5X87ncwaxZ85ztw6KEwYUJYdmfOnPI8v4jUmIKk8RU6KbGjI+wXzRMqV5BEAzSWLoXbbgsrR197Ley7L+y/P3z726qliCRa5qitOpW4IDGziWZ2s5n90sw+V+vyFBoky5eHZXImTAg/lzNIJk4M84/OPz/UUBYvhu9/PywYevnloflLRBIq1xIpdaaqQWJmPzazZWb2Utb2483sZTObZ2ZXdPcc7j7X3T8LnA5MqmR5C1FokETNWtE8oXKN3Gprg7337rxt9Gi45BL47/8OP2upe5EEU9NWF1OB4zM3mFlf4IfACcDewFlmtreZ7WdmD2dddkg95iTgKeAP1S1+V4UGSdTRXs4aydq1sHBh1yCJRCffevvt+K8lIjWSgCCpaunc/QkzG5u1+TBgnrsvADCzu4GPu/t1wIl5nudB4EEz+zVwZ659zOxC4EKA0aNHM2PGjJLKvH79+m4fu3nzHixdujMzZjzZ7fM89tiOwF68886LwH7MmvUSAwbE+4SfO7cJOISOjheZMWNFl/vffHMwcDhPPjmXwYOXxnqtSuvpOEt8OsbVUe7jfPCqVWx158U6/tvVQ8ztAizM+HkRcHi+nc1sMnAKMBCYnm8/d78FuAVg0qRJPnny5JIKN2PGDLp77DPPwD33wKGHTmbo0PzPM3NmuD7llP24+mrYbbd9KbFI73r99XB9xhn75VxaZ/XqcN3SMpHJkyfGe7EK6+k4S3w6xtVR9uM8eDDssENd/+3qIUgsx7a844zcfQYwo1KFKVbmXJLugmTJknCWzB12CD+Xo2mrrS2c8nf33XPfP2JE6HBX05ZIgnV0hH/kOlYPo7YWAZknF94VWFyjshSt0EmJb70FO+0UJqJC+YKktTV/86lZ6CdRkIgkmIKkIM8Ce5rZ7mY2ADgTeLAcT2xmU8zsljVr1pTj6XIqNEiWLAknLBswIFzKMWor14itbAoSkYRTkHRmZncBs4BWM1tkZue7eztwCfAIMBe4x93LMi/b3R9y9wtHjBhRjqfLqZgaSXTmy+HD49dINmwIfSQKEpEGl4AgqfaorbPybJ9ONx3n9ayYGslOO4XbTU3xg+Tll8OM9UKCZO7ceK8lIjWUgCCph6atRGtqCk1V3QXJunWhBhHVSJqa4jdttbWFa9VIRBrctm0KklqqRh+JWc+TEqPJiFGNpBxNW3PmhFUTxo/vfr/mZlixIrwXRSSBOjqgT31/VNd36WKqRh8JFB4kmTWSuEHS1hZmyfe0BE9zc3gfVjBLRaSS1LTVO/QUJNE6W+UOkp6atUDLpIgknoKkdyilaStOH8mmTbBggYJEpFdQkNRWNfpIoLAaSb9+MGpU+DlujeSVV0Kfh4JEpBdQZ3ttVbOPZN06eOed3PdHkxGj/rKmJli/vvQO8EJHbIGCRCTxVCPpHXqaSxIFSWT48HC9fn1pr9fWFt5XuRZqzKYgEUk4jdrqHXoKksxZ7RB/va22tjDsd+DAnvcdOjTspyARSSjVSHqHQmokUUc7lCdICmnWAi3cKJJ4CpLaqmZnO+QOkvZ2WLYsd9NWKSO3tmyBV18tPEggnM99RdfzXolIEihIaquane2QO0iWLw9rYpWrRvLqq+F9VUyQqEYikmAatdU7jBwZ/s65giR7VjvEC5JiRmxFFCQiCRUN7VSQNL4+fcKHda4gyZ7VDukgKaVpq60t9Hu0thb+mHoOEnf40Y9g/vxuTi8p0lt1dITrOh+1VfQy8mY2FNjs7h0VKE9i5ZuUmD2rHdJ9JKXWSMaNC6dxLlRzM6xcWZ9Nra++Cp//PAwZchBjxsCHPlTrEonUkShI6u0fN0uPMWdmfczsbDP7tZktA/4OvGVmc8zsm2ZWwGyGxldKjaTUICmmWSsqmzusWlX861XazJnhuqmpnY9+FO66q7blEakrCQmSQmokjwG/B64EXnL3bQBmNgo4Bvi6mT3g7ndUrpilMbMpwJTxPa21XgYtLfD88123L1kS+lAGDUpvGzQovC+Kbdpqbw8ntPrYx4p7XOakxOh2vZg5MxyfW2+dzQ03HMXZZ4fwveyyyr92R0cI83Xrwt9i3bqwmvJee4X5N5J87p1Po+CeumzZSt8332DEGy8ycJDBSSeFNuPItm3hsmZNGPZYKwnpIykkSD7k7luzN7r7SuA+4D4z62Ex89pw94eAhyZNmnRBpV+ru6atzNoIhPdrKeckmT8ftm4trUYC9dlPMnMmHHlkqJH89rdwzjlw+eWweDF84xvxm4bdYeFC+POf4dlnw/Xf/x6CY+PG3I8xC82H++6bvuyzT+iXGjAgXnkKKe/KlfDGG2F9tmHDwmXo0NCcmflZJ11t2wYvvghPPAGPPx6um5e30cY+vJ/HeZKj+SB/4EFOYgib3n3c1rHj6f/AvXDggfD734dvMi++GO7cYw/46U/hve+t/i/UKDWSKETM7HFgiruvNbPPAoOAm9x9S66g6W1aWkLT0datnc8R8tZbnftHIqUs3FjKiC2o3yBZvTqcoOvMM8PPgwbB3XfDv/0b3HBDOHa33174h3cUGs8/D3/9azo4li0L9w8YAAccAB/9KGy3XfgbDB/e+XrjxlCmOXPgpZfg4YfT/8sDB4bPmcMOS1/Gjy8+7NrbYdGisILz/Pnhknk737QnsxAoI0fCpElw1FHhcvDBPZ+XphjvvBPCNjoO3S1IOmgQjBiR+wKweXN4vs2bO1/MQlBGl75907ebmsLvOGJE+rpf6pPKPfyNVqwI7+cVK8LljTfgySfhqafC+wpgzBg4/nj43LJ74RF4gg/k/T36vz4PDjqIxZ+/hp1vurrznfPnw/veB//zP3DhheEPOG8eXHddCJkLLsj9T14OjRIkGUamQuQQ4ALgYeBW4NyKlCxhorkkK1Z0roEsWRI+cLKVcrrdKEj22qu4x9VrkDz9dLjO/KLXty98//uwyy7wla/A0qVw9tnh2/igQenL4MHhQ6WtLQTH88/DCy+kP0TMwnE64QQ49NDwN9h//8KWlfnEJ9K333knrLb84oshnP78Z/jxj0MZIf2hvvvuoUxDhnS+Hjw4fBC/9loIiwUL4B//CJ9Fkf79w+PHjQu1sz32CB+C27aFUzSvX5++bNgQnm/WLJg2LTx+8GA44gg4+mg4/HDYdVcYPTr83fN9/mzeHMJs4cJwmT8/HZ7z5qU/v/r1C8+TqybkHk5psHZtuF1JQ4eGy5o1+RdHbW2F006D978/HIsxY1J3nDm34NfpEiKZLroIjjsuhMnvf5/efv314VtLsd/wCtGAo7a2mlk/4FPA9e5+j5nNrlC5EidzUmIUJO5d19mKlNK01dYW/jmGDSvucVETb70FycyZ4f/jsMNgdsY7yQyuvDJ8ybvwQvjDH7p/nqFDQ0iccUaocRxwAOy3X3pQQxwDB4bn2m+/EGgQ/rfnzg2hEl3mzAnflDduDLXSbM3NISgOPTSUc9y4EB577BE++Ev5wvnWW/CnP4Vv4U8+Cddc03lF6WhY+ujR4bJhw77vBkh2LaNPn1C72mef8GEcNedNmNBzjdA9hNyaNZ0v0Dn8o8vAgeExHR0hUKNLR0dYuWHduvD41as7X69fH2onzc3hPR1dot8xOk1DF7NmhesLLoAbbwz9IX/8Y/EHHMIfLtvGjeFgVSJNG7BG8j3geUKT1hWpbUV+pDWuXLPb168P77F8TVvFrtxSyogtCN+OhwypzyA54ID8wXjeeXDKKeGDZNOmzk0jmzaF/7HW1vC/Xc0vbH37pvtOPvOZrve3t4fybdoU/v6jRqWHfJfTTjvBqaeGC4QP4OefD7XgpUu7XpYtG0Rrawiz97yn6yVzQEgxzML7uakphGJdcYf77gtVqwMPDNumTYPp0+Haa0M1FuD008Mf84gj6PjcxfS96+fFv9add6a/bZRLA3W2A+DuPzWz+4EOd99kZuOBWZUrWrLkCpJcs9ojTU2hSaFQHR2h3fqDHyytfPU2KbG9HZ55Bs7toWF0+PDKfAhXUtTOX44aUTGamkKfST4zZsxm8uTJVStPXTALbY+ZmppCtfCMM0Jb46hRoY0ype/UH8M+E+Hqq1kzZEembzyGs7iL2Wd+i0PGrsC+fl3u1/rkJ0MgxbVtW2i/Gzw4MTWSQuaRvNs66u7r3X1T6vY8d/909j69Va4gieaQ5KqRFNu09frr4Zt4qc2w9RYkL70Uamy1GAgj8q5x4zqFCBDa8q66Cv72N0YsamPi7J9x0sGLOPTuy/ngU//F8vOvyP1cEP5R47r00tCEEFW7IflBAjxmZl8ws90yN5rZADM71sx+Qp12uFdr9V9I90MUUyMpJkhKHbEVqbcgiSYiKkikbh1wAGy3HQce0pdpz+7CrbfC8239GTf1q/kf893vss9XvxpqQmZw773Fv+4tt4Tr229PTGd7IaU7HugA7jKzxWbWZmYLgFeBs4Ab3X1qBctYsmqt/gvpc7JnfljnmtUeiYKk0P65KEgmTiytfPUWJH/6U6ipvTuyRqSO9ekD//IvYULwWZ8ZzIl9pufe8Qc/oOXJJ9M/n346PPpocS82dmy4/tWvElMjKWQeyWbgJuCm1MTDZmCTu6+udOGSJntS4pIlYWhnrtEkw4eHptCNGwubRd3WFobElpqJ9RYkM2eG2ogaRSVJmptDheG1K09g9odOZ9KCe3p+0Ec+Uvg3xosvDp2hEALooYfC7ToPkkL6SC6Nbrv7Vnd/SyGSW64gGT06d6202PW2Sh2xFWluDqPEcg1NrbbFi0NTspq1JKl23x0mzf8FK789FYBH+3yE8Taf644psvaR7aabwvUpp4TrL385XCc9SIBLohtmdmbmHWY22sxOqNclUqotO0jyzWqH4paS37YtzFuIGyRQH2dKjIb1K0gk6UZdei64c8Di3/KJL47j1mf2L88Tn346nHxyYob/FhIku5lZNJDxR1n3/RQ4Ayhh0HXjyVUjydU/AsUtJb9wYZjRXGr/CNTX7PaZM8OktIMOqnVJRMpj9OjUBPc3RvPM5f/ddYfM5RLy2bIlXB99dJgVmrkMQwMEyUrgWjP7ONDPzN6fcd9O7n4e8JNKFC5pWlo6rzSab1Y7FNe0tWhRuI764EpRb0EyaVJhy5WIJMn228Omj72v6x3339/zg6PRpaefHtrDMxdQa4BRW6cBTxDW1zoV+L6ZfcrMvgQsA3D3X1euiMnR0hIGWaxaFSbcLV9enqatqJazww6ll61egmTzZnjuubAGnkhDKnUESbRQXDSvJTNI6rxGUsjM9q3AL939XgAzex24iLBUSsWXZ0+SzEmJW7aEgRrlaNqKVq9thCB57rnQ4a/+EZEsUY0kGpqZoCAppEZyLvCcmd1tZucBa9z9Mnf/vLu/VtniJUtmkHQ3qx2Ka9qKgiTOSanqZeHGaCLikUfWthwiFfXP/1z8YxJcI+kxSNz9s+5+MPCfwHbAVDObZWbXmtn7zaxuf8NqzmyHzkHS3ax2KL5pa8SIeH0KAwaEWlA9BMn48fFqVyJ172c/67qtp7kkURt2VCPJXHY56UEScfe/u/uN7n48cCzwFKH/5JlKFS6uas5sh9w1knxBMnRoaEottEZSjg/eWk9KdE9PRBTpdbqbxPX22+mVgxPYtFXMMvIAmNlQYLO7TwfyrBPQO0VNT8uXp/vb8gVJtPR2bwqSBQvC76IgkV5py5b8J3fJnDeQq2kr6aO2zKyPmZ1tZr82s2XAy8ASM5tjZt80sz0rX8xkGDgwNB9FTVsjR3Z/jodCz5K4bFm6thNHc3NlJyRu29Z97V0LNUqvlu/Ujtn3Re3eCaqRFLT6L7AHcCWwo7vv6u4twNHA08DXzayEnqXGFE1K7G5We6TQpeSXL09GjWTKFDj22PyvMXNm+J0rcUZSkbpz+eWdf47W0Mpl8+Zwfd996dpHgoKkkKatD7l7l8Y9d18J3AfcpyVS0qIg2bQpf7NWpJCmrW3bkhMkf/lLqIkdcQQ8/HDXc8vPnBnuq/P/CZHy+Na3wofAF78Yfj7qKHjxxXBqzWxRkGSu8JqgIClk1NZWADN73MyGp25/1sz+j5kNyNxHOtdIyhEkK1eGMClXkKxfn37PlpN7CKkTTgi/05FHdj7X+tq14X9IzVrSq0TLXETefDP3ftE/ZebQzEYKkgwj3X2tmR1CmIi4HXBrZYqVXFGQLFnSc9NWIX0k0RyScvWRQGX6SVavDrP5jzsunEJ3113h+OPh1tQ75JlnQtgoSKRXyQ6SfJ3tUZBkdqomaPhvMaO2tppZP+BTwPXufo+Zza5QuRKrpSXURrqb1R4ppI+kHMujRDJnt++yS/zny5QZeGPHhhNXnXEGXHhhOBnQkCFhpNrhh5f3dUXqWqFBEnW2ZwZJgkZtFRMk3wOeJyyNEp20eFjZS5RwLS3pkUuF1Eh6CpJyLI8SqeTs9uzAGz48nJPnssvghhvCGST32y+9NIxIr5AdJP3yfOTmqpE0YtOWu/8UOBzY1903mdl4YFbFSpZQmU1QhfSRrF3b/ZDZSjRtVTJIMsvZrx9873vwgx+E/6djjin/64rUtewgyf450uhBYpZeytLd17v7ptTtee7+6ex9ertigmT48NCv0N3w8mhyY1SbiKPaQRK5+GJ49VX42tfK/7oide3cczv/3N6ee79GDxLgMTP7gpntlrnRzAaY2bFm9hPCwo5C5w/SQpq2oPvmrWXLQojkqxEXIxpZWIkg6anmNG5cYeemF2koY8aEkSaRXhwkxwMdwF1mttjM2sxsAfAqcBZwo7tPrWAZS1btRRsh/UHav3/nIeG5FLJwY7mWR4EQRtttV7kayfDhOlmVSBeZ/xT5giRqlsjct5FGbbn7ZuAm4KbUxMNmYJO7r6504eJy94eAhyZNmlS186ZEQbLjjj2f36aQc5KUa3mUSKUmJS5fXt5yijSMzEDo6Mi9z+bN4dtn5uisBI3aKqSP5NLotrtvdfe3khAitTJkSLj01D8ChTVtlWtWe0RBIlJlmUHSXdNW9sJ8Dda0dUl0w8zOzLzDzEab2QlaIqWzHXYoLkiq1bQFlQuSctecRBpGriBZuTIs9RDpBUGym5mlPvL4UdZ9PwXOAH5e1lIl3He/C1df3fN+PTVtbd0a3m9JCJJy15xEGkZmv0fUtHX44bD//untPQVJOUbbVFAhpVsJXGtmvwf6mdn73f2J1H07u/tHzOxjlSti8px0UmH79dS0FX3gV6KPxL3nPpxCRetsqUYikkOuGsm8eZ33eeed7oOkzmfyFhIkpwGjCetrnQp838xuAHYElgC4+68rVsIG1lPTVjmXR4k0N4cvPxs3lm847po1ofYWgwVWAAAQHElEQVSkIBHJobs+EneYPh3uvDM90SuSGSR1PlWvkFFbUe3jXgAzex24EDgQqMA6sr3HsNQCM/lqJOVcHiWSOSmxXEHS3WREkV6vuyDp6IB77w23s9uc863LVYeKGlNmZgcC5xNqJsOAoypRqN6ib9/wYV6rICmXSpRTpGFk9m9kD//dsgUOPjj34/onZwxTIcN/J5jZV83s78BtwArgA+5+OKH/RGLobin5cq6zFalEkKhGItKD6HwK2TWSLVvSQfPhD3e+L0FBUkgfyd+BZ4FT3f2lrPu6WW5QCtHdUvLLl4f32MiR5Xs9BYlIDUyZEq5zBUm0kOPPswa/RkESdabWsUKatj4BvA78zsx+llp2JDlRWee6W0o+mptRzkmtChKRGohqHbmatqIgyf5Hjx6TgGWzC+lsfwB4wMyGAicDFwG3mdl0oL7HpCVAT01b5f5wHjkyvF/L3UcybFjX0YsikhKFQnc1kuwgaW6GJ57I34dSR4o5H8kGd/+5u58ITASeBl7s4WHSg56atsrdgd2nT1hNuNw1EnW0i3SjlCABOProRCybXVKjibuvdPf/cff6r3PVuZ6atirxAV3u2e1aZ0ukB9ESJ8UGSUIkt+QNoqemLQWJSAPI7CPJDBMFiZRDvqatzZvD9kp8QJc7SLRgo0gPMmskGzemtytIpByamsIyO1u2dN5eieVRIuUMEnf1kYj0yCwERXs7bNiQ3v6+90FbW7itIJFS5Vu4sZKzxZubYcWKEAJxrV2rdbZECtKvX9caCcDjj4drBYmUKt9S8pUOkq1buz+hVqE0h0SkQP36hT6SzBoJpNfUUpBUl5kNNbPnzOzEWpclrp5qJJXqI4HyNG8pSEQKFNVIsv/ZoyCp8xV+u1PVIDGzH5vZMjN7KWv78Wb2spnNM7MrCniqLwP3VKaU1ZVvKflK95FAeYJECzaKFKhv3xAk2f/sAweGEElwkFT7tFtTgR8QzqwIgJn1BX4IHAcsAp41sweBvsB1WY//DLA/0AY0xDzq7pq2Bg1KLzVfTqqRiNRA1LS1Zk3n7X36JLpZC6ocJO7+hJmNzdp8GDDP3RcAmNndwMfd/TqgS9OVmR0DDAX2BjaZ2XR335ZjvwsJ501h9OjRzJgxo6Qyr1+/vuTHFuK114YCh/L003MYNGj5u9tfeGEvhg8fyeOPP13211y8eBBwBE8+OZchQ5bGeq5nntkNGMfcuU+wYEGXP0PBKn2cRce4WvId5yO3baPj4YdZ3KcP4zP3X7WKIWY8keS/jbtX9QKMBV7K+PlU4LaMn88BflDA85wHnFjIax5yyCFeqscee6zkxxbi9dfdwf222zpv/+hH3WMUu1tr1oTXvOGG+M916aXuQ4fGf55KH2fRMa6WvMc5DJTsepk40X3gwKqWsVDAbC/gM7Yeziifq2Gwx4Gp7j61/EWpvu462yvV79DUFFaoLlcfiZq1RGLYsiXxTVv1UPpFwHsyft4VWFyjslRdLYLErHyTEjUZUaQEmSet2rpVQVIGzwJ7mtnuZjYAOBN4sBxPnDp3yi1rsju36kj//qFTPXMgh3vlv+mXM0hUIxEpUuY/jWokxTGzu4BZQKuZLTKz8929HbgEeASYC9zj7nPK8Xru/pC7XzhixIhyPF3FZK8AvGFDWGurkt/08wXJpk1w//3w6KOFPY+CRKQEo0albzdAjaTao7bOyrN9OjC9mmWpJ9kLN1ZjbkZzM7yUms3zzjvwyCNw993w4IMhyEaOhJUrux/aXo2ak0hDyvxym3ne9oRKdukbRPZS8tUIku23h0WL4LzzYNq0MLR91Cg4++wQELfdFsoxenT+51i3LvwPqI9EpEjDM04uu2VLenZ7QiW7PtWDJPSRQNemrUoujxLZddfwmtOmwT/9E/zmN7BkCdxyC5x6atjnlVe6fw5NRhQp0ZAh6dsN0LSV7NL3ICl9JNlNW5VcHiXyhS/AjBmwdCncfjscf3x6IMmECeH65Ze7fw4FiUiJ+veHhQvD7W3bFCQSXy1qJMOHwwc+EJb5ybbbbmG7gkSkQvr371wrUZBIXLn6SJqaYPDg2pSnb1/Yc8+em7a0YKNIifr3T581ERQk9SwpfSS5Rm3V+lv+hAmqkYhUTP/+nUdqKUjqV1L6SJqawknTOjrCz/UwW7y1FebPD/2A+SxfHmrnmTV0EcnjnowzXyhIpNyyl0mp5PIohWptDadOeO21/PtoMqJIEU47DcaMCbfVtCXlln1Oknpo2mptDdfd9ZPUQ+CJJEpUC1GQSLll1kjc66Npq5AhwKqRiBQpM0jM0gGiIJG4Mk+3u3p1aFKqdZCMGhWWUVGQiJRRFBhRoETXCpL6laRRWxBqJPU0pLa1NX+QRDUnBYlIETx1qqVo9q+CpP4ladQWdA6SeviAbm3N30eyfn1YobgeyimSOFGQRP0kChKJK7NpqxrLoxRqwoSw/lbmZMlIPZVTJDFUI5FKqeemLcjdvKXJiCIlyA4S1UikXHI1bTU31648EQWJSJk1aI1E5yOpAwMHhvfV2rWh72G77Tqf0rlW9tgjvL9z9ZMoSERiUI0kOZIyagvS623VwxySyIABsPvuuWsk9dQEJ5IYDVojSXbpe5CUUVuQXkq+3maL5xsCvHx5WJ146NDql0kksRQkUknRUvL1sDxKpmgI8LZtnbdrDolICdTZLpVUj01bEIYAb9oEb77ZebuCRKQEqpFIJTU1heVR3n67voIk38itemuCE0kE1UikkpqawpLt7vX1AZ0vSFQjESmBaiRSSU1NsGpVuF1PH9A77QTDhnUdAqwgESlBFCRatDE5kjb8N1JPNRKzrqfd3bAh9JsoSERKpKat5Eja8N9IPQUJdB0CXE8LS4okipq2pJLqPUj+8Y9QCwEt2ChSMnW2SyVFTVt9+oSTStWT1tbw/p8/P/ys5VFESqQaiVRSVCNpbq6/91T2aXcVJCIlUpBIJUVBUo/NRQoSkTJR05ZUUtS0VY9BMmwY7LJLOkiWLYNBg8J2ESmCaiRSSVGNpF6/5U+YkJ5LEs0hMattmUQSJwqSqCaiGomUUz03bUF6CLC7JiOKxBZ9C4tqJAn/VqYgqRP13LQFIUhWrQprgSlIREo0YEC4joIjauJSjaR+JWlme0sLXHQRnHRSrUuSW+aaW1qwUaREjz4KV10Fo0eHnwcODNcKkvqVpJntffrAzTfD/vvXuiS5RSO3XnlFNRKRku21F1xzTbpGEtVQFCTSG4wdG97zf/0rbNyoIBEpCwWJ9CZ9+8L48fDUU+FnBYlIGahpS3qb1lZ44YVwW30kImWgGon0NhMmpM/drhqJSBlENZJofklCKUikYNHILVCQiJRFVCNpb69tOWJSkEjBFCQiZRbVSBQk0ltEQTJgQOfzp4hIiaIaydattS1HTP1qXQBJju23D+dKGTIk8Ss6iNSHBmnaUpBIUfbaC955p9alEGkQDdK0pSCRotx4I2zZUutSiDQINW1Jb3TYYbUugUgDaZAaiTrbRURqpUFqJA0dJEla/VdEeqEG6Wxv6CBJ0uq/ItILqWlLRERiUdOWiIjEohqJiIjEohqJiIjEohqJiIjEohqJiIjE0r9/uFaNRERESqIgERGRWPqlVqlS05aIiJRk+HB43/vgzjtrXZJYtGijiEit9OkDTz1V61LEphqJiIjEoiAREZFYFCQiIhKLgkRERGJRkIiISCwKEhERiUVBIiIisShIREQkFnP3Wpeh4sxsObAayD55+4gCtjUDb1eudD2Wp1KPL2Tf7vbJd18hxzTXtqQc52If29P+5TrGubYn9RgX+3i9lyv3+DHu3tLjs7h7r7gAt5SyDZhdyzJW6vGF7NvdPvnua/TjXOxje9q/XMc4zzFN5DEu9vF6L9fu8dGlNzVtPRRjW7XEfe1iHl/Ivt3tk+++Rj/OxT62p/3LdYxzbU/qMS728Xov1+7xQC9p2orDzGa7+6Ral6PR6ThXno5xdfTG49ybaiSluqXWBegldJwrT8e4OnrdcVaNREREYlGNREREYlGQiIhILAoSERGJRUESg5mNM7P/NbNf1rosjcTMhprZT8zsVjP7ZK3L06j0/q0OMzs59V7+lZl9uNblqYReGyRm9mMzW2ZmL2VtP97MXjazeWZ2RXfP4e4L3P38ypa0MRR5vE8BfunuFwAnVb2wCVbMcdb7t3RFHudpqffyecAZNShuxfXaIAGmAsdnbjCzvsAPgROAvYGzzGxvM9vPzB7OuuxQ/SIn2lQKPN7ArsDC1G4dVSxjI5hK4cdZSjeV4o/z1an7G06/WhegVtz9CTMbm7X5MGCeuy8AMLO7gY+7+3XAidUtYWMp5ngDiwhh8jd695edohV5nNuqW7rGUcxxNrO5wNeB37j7X6pa0CrRP2lnu5D+JgzhA22XfDub2fZmdjNwkJldWenCNaB8x/t+4BNm9iNqu/xEo8h5nPX+Lbt87+cvAB8CTjWzz9aiYJXWa2skeViObXlnbLr7CqAh3xhVkvN4u/sG4NPVLkwDy3ec9f4tr3zH+XvA96pdmGpSjaSzRcB7Mn7eFVhco7L0Bjre1aHjXB299jgrSDp7FtjTzHY3swHAmcCDNS5TI9Pxrg4d5+rotce51waJmd0FzAJazWyRmZ3v7u3AJcAjwFzgHnefU8tyNgod7+rQca4OHefOtGijiIjE0mtrJCIiUh4KEhERiUVBIiIisShIREQkFgWJiIjEoiAREZFYFCQiIhKLgkRERGJRkIhUiJkNNrPHU+epwMzczH6WcX8/M1tuZg+X8NwDzOwJM9PCq1JzChKRyvkMcL+7Ryfn2gDsa2aDUz8fB7xZyhO7+xbgDzToGfckWRQkIiUys3PN7Dkze8HMnsyxyyeBX2Vt+w3wsdTts4C7Us811sz+njpX/Qtm9kszG5LxWp9KbX8+o1YzLfUaIjWlIBEpgZk1AV8GjnT3/YEpWfcPAMa5++tZD70bONPMBgH7A89k3NcK3JJ6vrXA51PPtQ9wFXCsux8A/Ftq/5eAQ8v5e4mUQkEiUpoOYDBwg5lNcvfVWfc3A9nbcPcXgLGE2sj0rLsXuvufUrfvAI5K3T4W+KW7v516jpWp6w5gSyrURGpGQSJSAnffCOwL/Am4xcw+n7XLJmBQnoc/CHyLVLNW5tPm+dly3BcZCGwupMwilaIgESmBme3p7hvc/W7gYbJCw91XAX1TTVjZfgz8P3d/MWv7bmZ2ZOr2WcBTqdt/AE43s+1Trz0qdb09sNzdt5bllxIpkYJEpDRXmdnLZvYXYHfgphz7PEq6eepd7r7I3b+bY/+5wLlm9gIwCvhRav85wNeAx83seeDbqf2PoWvzmEjV6cRWIhViZgcBl7n7OQXsOxZ42N33LeL57weudPeXSy6kSBmoRiJSIe7+V+CxaEJiOaVGhU1TiEg9UI1ERERiUY1ERERiUZCIiEgsChIREYlFQSIiIrEoSEREJBYFiYiIxKIgERGRWBQkIiISy/8HGtMK5FYYFU0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEeCAYAAACt7uMeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztnXmYHFXV/z8nM1km+wwJIWRhAmELkYAkQJBXAgJhX5QtgoCyqICigg/rw4s7/AR5ZRMBJSAIBgQkLBIBAxISZYeEEIghZIXse8gyc35/3C66pqeX6u6q7uru83mefqq71jN3uutb55x7zxVVxTAMwzCKpUO5DTAMwzCqAxMUwzAMIxRMUAzDMIxQMEExDMMwQsEExTAMwwgFExTDMAwjFExQjIpGRMaLyM9DPufZIvJymOfM49rNIqIiUp/4/IyInOXb/nMRWSYinyQ+nygi80VknYjsXQ6bDcPDBMWoCERksoisFJHO5bbFT9Tio6pHquq9iWsNAi4BhqnqdoldbgAuUtXuqvpmVHYYRhBMUIzYIyLNwP8AChxXVmPKyw7AclVdkrJuRiEn87wgwwgLExSjEjgTmAaMB85Ks72PiPxDRNaKyIsisgOAOG4SkSUislpE3hGR4YltvUTkPhFZKiIfi8jVItLu95Aagkqsmywi54rI7sAdwOhEyGlVYntnEblBROaJyKcicoeINKT7w0SkLrHvMhGZAxydst271qHAP4DtE9d6UETWAXXA2yLy38T+24vIXxN/10ci8n3fua4VkUdE5H4RWQOcLSIdRORyEfmviCwXkQki0pTyt5+V+FuWichVKbZfmTh2rYi8nvCiEJHdEv+TFSIyS0ROyfofNqoCExSjEjgTeCDxGisi/VK2nw78DOgDvJXYD+Bw4MvALkBv4FRgeWLbLUAvYEfgoMQ1vpmPUao6E/gOMDURcuqd2HR94pp7AUOBAcA1GU5zHnAMsDcwEjgpw7WeA44EFiWuNU5Vuyc2j1DVnRKCOBF4O3HNrwA/EJGxvlMdDzyCa48HgO8DJyTaYHtgJXBbyuUPBHZNnO+ahJAC/AgYBxwF9AS+BWwQkW448fszsG1in9tFZI8MbWBUCSYoRqwRkQNxYZ0Jqvo68F/g6ym7PaWqL6nqJuAqnMcwCNgC9AB2A0RVZ6rqYhGpw4nLFaq6VlXnAjcC3wjBXsGJxA9VdYWqrgV+CZyW4ZBTgP9T1fmqugL4VRGXHwX0VdWfqupmVZ0D3JVy7amq+riqtqrqRuDbwFWquiDRftcCJ6WEw36iqhtV9W2cWI1IrD8XuFpVZ6njbVVdjhPIuap6j6puVdU3gL+SQSyN6sFiqEbcOQuYpKrLEp//nFh3k2+f+d4bVV0nIiuA7VX1BRG5FffEPVhEHgMuBRqATsDHvnN8jHuqL5a+QFfgdactAAguNJWO7f32p9iULzvgQmKrfOvqgH/5Ps9vewg7AI+JSKtvXQvg9wI/8b3fAHie0SCcwKezY78UO+qBP+X8C4yKxgTFiC2JvMMpQJ3XTRboDPQWkRGJJ2ZwNzbvmO5AE7AIQFVvBm4WkW2BCcCPcU/hW3A3vvcShw4GFqYxY31i2RVYk3i/nW97arnuZcBGYA9VTXe+VBb77U/YUSjzgY9Udecs+6TaOx/4lqpOSd0x0Rki1/V2AqanWf+iqh6W43ijyrCQlxFnTsA9LQ/D5SP2AnbHPXGf6dvvKBE5UEQ64XIp/1bV+SIySkT2E5GOOGH4DGhR1RacuPxCRHokkvg/Au5PNUBVl+KE5oxEEvpbuJuox6fAwMS1UdVWXJjppoSIISIDUvIYfiYA3xeRgSLSCFyedysl+Q+wRkQuE5GGhL3DRWRUlmPuwLWD15Ghr4gcH/B6dwM/E5GdEx0g9hSRbYAngV1E5Bsi0jHxGuXLvRhVigmKEWfOAu5R1Xmq+on3Am4FTvfF+f8M/C+wAtgHl6QHlyi+C5do/hiXkL8hse17OJGZA7ycOMcfM9hxHs6zWQ7sAbzi2/YCrtvuJyLiheUuA2YD0xK9qZ7DJbXTcRfwLC438QbwaLYGyUZCKI/FCe9HOG/pblzng0z8FngCmCQia3G96fYLeMnf4ARxEs57+wPQkMgbHY7L3SzChcyux3mXRhUjNsGWYRiGEQbmoRiGYRihYIJiGIZhhIIJimEYhhEKJiiGYRhGKJigGIZhGKFQUwMb+/Tpo83NzQUdu379erp16xauQUYbrI1Lg7VzaaiWdn799deXqWrfIPvWlKA0Nzfz2muvFXTs5MmTGTNmTLgGGW2wNi4N1s6loVraWUQClwOykJdhGIYRCiYohmEYRiiYoBiGYRihUBOCIiLHisidq1evLrcphmEYVUtNCIqqTlTV83v1ylYjzzAMwyiGmhAUwzAMI3pMUGqA2bNhfuo8fYZhGCFTU+NQaonVq+Evf4Hx42HqVNhvP5g2rdxWGbVGa6t7mPngg/avBQugSxfo0QO6d2+77NsXDjoIDj0Utt++3H+FERQTlCqipQWef96JyGOPwWefwbBhsMcesDDIZLSGUSDr1sGsWe71/vvJ14cfuu+hR7dusMsuMGoUfO1rsGmTO3btWrdct84J0CuvwB8T050NGwaHHebE5aCDnOAY8cQEpUr4+9/h3HOdcDQ2wjnnwNlnwz77wCWXwF13ldtCo5pobXUe78MPw9/+Bh99lNzWoQPsuCPsthscfjjsuivsvLMTkv79QSTY+d95B/7xD3juOfj97+G3v4X6ehg9GsaOhSOOgL33dtcz4oEJSpXwwAOwYQNMmADHHQedfZOt9u7tnvy2bnU/SMMoBE9EJkyARx5xDy+dOrmb+7nnOgHZbTfYaae2379C6NAB9trLvX78Y+flvPKKE5dnn4Wrr3avvn2daI0d65b9+oXztxqFYbeXKmHePBg+HE4+uf223r3dctUq6NOntHYZlc+aNfDLX8L99ydF5Igj4Prr4dhjoWfP6G3o0gUOOcS9fvlLWLIEJk1y4vLss+6BCuCLX3S2HXGE82TsAaq0mLNYJcyfD4MGpd/W2OiWq1aVzh6jOnj8cdh9d/h//8+FT++/H5YudWGu008vjZikY9tt4Ywz4E9/gk8+gddeg5//3OVorr8evvxl9/B00knwhz+4DgBG9NSEoFT7SPnWVveDySQofg/FMIKwcCF89atw4onuxjxtWvlFJBMdOjixu+oqeOklWLbMheROPtnZfe657rdx7LGu92NNsWGDixl++mlJLlcTglLtI+U//RS2bIHBg9NvN0ExgtLaCrff7rySZ55xT/uvvQb77ltuy4LTu7frQXbXXc5zf/dduPZa2P7pu5m/3SjmfqTlNrF0PPQQ3HADXHppSS5nEcYqwBu0mMtDWbmyNPYYlcmMGXDeeW7c0qGHwh13uAR7JSPicovDhwPXngefwV6jFnLHkwPZf/9yW1cCnnzSLZ97riS9cmrCQ4mK1la4997Cxnhs2hSeHfPmuWUmD8VyKEY2NmyAK65wPao++ADuu88lvCtdTDLx1vJBnHXQXCZMKLclJeDNN93I0KuucoISMSYoBaIKF13kxnp4A7CCMnOmGxH83nvh2BLUQ4mzoLz6Kpx55r68+mq5Laktnn7aDXy97jqX5H7/ffjGN4KNFakoFi9u8/Hhzmfw4qm3cefFM9At0d9oy4IqLFrk/rEXXeS6ykWMCUoBqML3vge/+537vHx5fsd/8IF7WPAPBiuGefOga9ekJ5JKt25QVxfvkNdzz8H8+V058sjwhNbIzMKFLml99NHQ0ACTJ8M991Rxt/KNG+HUUz//uOfaKdzGRZx/83A+3m4/tkwtbGrwWDJjhusz/dJLsHlzSWvXmKDkiSr84Adw221uBPrgwfnfqFescMv168Oxaf58Z0emp0oRJzZx9lCmT4devTZTX+8GqM2dW26LqpOWFnjkkQHstpsLr//iF/DWW66kSVWz444uQd3a6kZJHnjg55uaV7xBxwNGhRuHLhevvuoSRtOmgTefvQlKPFGFH/4Qbr7ZLX/9a9hmm6RABMUToA0bwrEr2xgUj9694y8ou+22lkmTnNAedpgbX2CEwyefuLEkw4bBbbftzJe+5Nr8yivdQMWaQcQN4588uf22EoSEIidddzwTlPih6jyS3/4WLr4Ybrwx+eRfbg9l3rzMCXmP3r3jG/LautXF7ocMWc+ee8JTT7nQ79ix8RbBuLN1q/NCTjgBBg6Eyy5zpUquvXY6zzxTvUn3QNTVpV2t60N6yosT/fuX7FImKAFQhd//fkduusnltm66KRleamrK30MJU1A2bXJPn7k8lDiHvGbPdqHe5mbXIAccAI8+6jovHH10eMJbC2zYAP/+t/M8Bg92g/mmToUf/ci158svw0EHLau+pHshvNY+b7Jm932rb/TjgAElu5SNQ8mBKlx+OfzlL4O54AIX7vL/GBsbCw95hXGj9LosBwl5xXWSrenT3XLIkGSDjB0Lf/6zy6OedJIbpV1ToZkALF3q8h9vveV6h771lisf39rqRo8fdZSrOn300dCxY7mtjSH77OO+XI888vmqXvNnsOXhx+l47lllNKxAbrml/bpFi4qv1JkHJig5WLsWJk6E445byC23DGj3ZNfU5ARCNXhXyzA9FE8kKjnkNWOGa7sddmgbbjjpJFe2/Lzz3JP22LGw3Xbu1a+fWzY1uWNbWpwHtmKF+zu95fr1rorA1q1tX1u2uONHj3Z5hQwRkNixfDk8+KCb8+b115PrBw9240hOPtktR492f5+Rg/HjXUPdeuvnq17+23IOPrd8JhXM97/ffl2Ju+3VhKCIyLHAsUOHDs372J49YcoUePPND+nQob3r2NjowjUbNrjuuUGIQlAqOSk/fbpX8ry13bZzz3XtdPnlbrBdKh07um6va9YUfv0ePVwuc/Ro99p/fydUcWHrVldR95574IknnBjuvbcrizJypBOQONlbUXTr5r5kPkE5+MlLWH3fTvQ68/gyGpYnmqGcTIld05oQFFWdCEwcOXLkeYUc39iYeRIf74e8cmVwQQmzl5c3Sj5IDmXTJtdjMm6dWaZPT5TGyMDFF7uHr1WrXN2yTz5pu1y/3v19TU3tl926ud9Ufb17ee/r6lzX5KlT3WvaNPjVr5ynA3DwwS4P8ZWvlG+Q38cfu+7pXkXdPn3gwgvdYNoRI8pjU1UyYoT7MR5wgIsbAr3OOgEmne6+FLl+XHEgrC6jRVITghIlnqCsWOF60gQhbA9lm23cwMZs+EfLxykUsmmTmyb2pJOy7+f1qGtsdJM4hcHQoe71jW+4z+vWuTztSy+5UNthhznP5corXcitVDMDzp/v5vz4wx/cg+fRRzsROeooyyNFRkODS0Q9+CCzfnw3uy58wU2y8sADbj7tE04ot4XZSRfP3n33kpthvbyKxBudHjQx39qaDD2FISjz5gV7gIprgchZs5xXsMce5bbElcMZMwauuQbmzHGismyZu5eMGOE6CURZDmnhQteLcOhQJybnnuvsePxxZ4OJSQkYN45tpj3Vdt2JJ5bHlqCoOk/Kz4UXlqXkhAlKkfhDXkFYvToZ7gzLQ8mVkIf4Foj0enhlC3mVg86d4fzzneDdf7/7n51+upsf/dZbi8vZpLJ4sQvr7bSTE7Gzz3ZdqW+/vTKiLdVGn4FdmLfTmLYru3cviy2BeOUV92Xx+NnP0vf4KgEmKEWSr4fi368cHkrcBGXGDJfT2GWXcluSnvp6JyTvvOM8hb59XR23AQOcNzFzZmHn3bzZDTr8+tddVZDbbnM1/D780IlKkIcEIzr6zfgnrfiSZ+vXZ058l5t169p+bmkpW+LPBKVI8vVQPEFpaCg+j7Z2rfN4gtx84hrymj7diUncwzkdOsDxx7vk/b//7WYzvOsu1+X40EOd2OQKh7W2wr/+Bd/5jhu8fOyxrvfWN7/pPKG774bm5pL8OUYOOneGDqQIyGmnuSeBuOGfjXHPPV24q0xYUr5Iund3T7FBPRTvhj5oUPEeStAuwxDvkNfIkeW2Ij/23de9brjBicDvfufC7AMHOnHs1s11kvC/Nm1yud35893n44933snhh8dfTGsVPeUUxD9pyoQJ7mng4IPLZ1Q6/KX5//nPsvYhN0EpEq/3Ub4hr0GDXKeSYsg1sZYfb/bjOAnK+vWuhP9ZFTgoGVz464or3JTdEye6iamWLHH/4w0b2r62bnW9xn71KycmcQ7JGw657z7eW7Edw567Obly48byGZQJv6D07Fk+OzBBCQVvtHwQPEEZONDVVSqGfDyULl3cK04hr5kzXVg6bgn5fKmvdx5Kts5A+VRSMGJC584MHd0HnvOte+opNziphOVMsvLKK65iLbgKyhFP8ZsLy6GEQD4eindDHzDAhUG8gXSFMG+ei+0HrU4dt9Hyce3hFQUmJpVJp9O+1nbF7bfDpZeWx5h0eNPFduoUi0ltTFBCIF8PpVu3ZE6jmDzK/PlOTII+lMSt4vCMGe5Br6bLqBvxZtgwZn+ojOPPyXVTp5bPnlQ+/NAt083vUgZMUEIgnxL2K1YkS4JAcT29gkys5SeOHsruu1dOYUajNhk6FFYdMS65IldZilLR2gpvvOH6r48eXW5rABOUUMhnkq2VK93+nqAU46EEmVjLT9wqDueq4WUYceGCC3wfFi6MtmRCUD791I1BKUOJlUyYoIRAU5N78g+SD0n1UAoVFNX8PZQ4hbxWr4YFC+JRcsUwcnHUUXD2tk+zsMtOrh7OKaeU16A33kgmT0s4xW8uTFBCIJ8xHitXhiMoS5e6pH6lhrxmzHBL81CMSqCuDnb74ZEc/VliMq7HHks742PJeOCB5HsTlOoin9HyK1Y4AfLCsIUKStCJtfx4ghKHChK11MPLqA7OOQdmdtqLK775iVsxalTbG3sp8c8TX8I543NhghIC+dTzCivkFXQeFD+NjS4sl1r6pxzMmOHawGpWGZVC374u0nXbw9smV55xRnncfv84mBjNR1ETgiIix4rInatXr47k/EE9lI0b3QRXYfTyKtRDgXiEvaZPd/mTUs0xYhhhcOGFsHZdyqAi74mylPifCks8K2M2auLnrKoTVfX8Xl79kZDxT7KVDU9wwujlNX++G/mez5TRcSoQOWOGJeSNymO//eCLX4QPO5f5y+vdOGL2I6oJQYmaoCEvb3tYIa+BA/MbgR0XD2XpUtfj0fInRqUh4roQH7jpeTYM3Dm5YcqUcC+UK9G5bp0r0OclI2OCCUoIeIKS68nf2x6GoASdWMtPXCoOWw8vo5IZNw429+7HC12PSa488EA3K1oYfPyxiwX/5S+Z91m/PnkTiREmKCHQqZP73wb1UBob3TEdOhTnoeQ7m19cPBRPUGLmrRtGILp2dXPYfPBhyu0zrN4us2a55e9/n3778uVuUp4Ylqw2QQmJIPW8/CEvESdChSTlt251Favz9VDikkOZPt3ZEqPu84aRF9/9LtTplrYrP/ssnJN7tYgyPfkdfDC8+655KNVMkHpe/pAXuO9DIR7KokWujE++Hkpc5kTxSq5YBV6jUtl5Z9javHPblWvWhHPytWvdMlOv1HffdcstW9JvLyMmKCERpIT9ihXu4aNHD/e5UEEpZAwKuKrEPXqUV1BUrYeXUR3Ufe8CfsI1yRVjx4Yz0NETpjlz4De/ab/d80yCVqQtISYoIRE05NXYmHwyL1RQChmD4lHuApGLF7vrW0LeqHROPrUDj5Myq9q99xZ/Ys9DAbjkkvbbt00MrDRBqV6CeCheHS+Prl1L66FA9AUit26Fiy+Gf/wj/Xavl6N5KEalM2AA7DIqZWxbGDM5+kNnX/pS++39+rllWDmbEDFBCYl8PBSPYjyU3r2TobN8iLpA5OzZcPPNzvu/5pr2FZity7BRTYw9JUVQwhi1vmaNi08fd1z6vIx3E4lhmYn4WVShNDW5B4aNGzPv49Xx8ii0l1e+Zev9RB3y8sTqi1+En/0MDjsMPvkkuX36dOex9+0bnQ2GUSqOOaM3UzgguaJTp+JPunYt9OzpfiTLl7ff7iXj33ij+GuFjAlKSAQZLe9NruVRTFK+0KKKUXso3rlvuQXGj4dp02CvveCFF9x6r4aXYVQD227XgZ8c5hsl/9e/FnejX73aDWjs2RO22QaWLWs/an7TJhgzBkaMKPw6EWGCEhJBCkSm81AKDXkV6qFEnUPxzt27N5x1Frz6qvubDzsMfvpTeO89C3cZ1cVpp/k+bN0K++xT+MmuvtqJyNy5rlDf5s3tbxKbN4fjCUWACUpI5PJQWlvdzbZYQdmwwXnBxYS81qwJNrtkIfgFBZw38p//wOmnw//+rxtMbIJiVBMnnph7n8B4gxohWfl12bK2+2zaFE7yPwJMUEIil4eyerXzXP0hL6+XVz4TXhXTZRiSN/qIKvl/Lij+v7N7d9eb8u67Ybfd4JBDorm2YZSDxkZ4cocLwzmZ14OrXz8X8gITlFokVwl7f9kVj27dnKeweXPw6xTTZRiiLxC5cqX7rnfp0na9SGLGu5kwdGg01zaMcrHq57eymBAmuvJuBu+9l7xZpP5YTVCqn1wVhzMJCuTX0yssDyUqQVm1KnkNw6gVjjsOWqhPrih0jMjmzS7s1dTkEvPQvuuw5VCqn5493fcgk4fin1zLo5AS9vPnu6f9AQMKszPqApEmKEYt0rMn3PCVZ5hdt4tb8emnhZ3ILxbpBGXKFFiwwDyUakfE3UjzDXlBfoIyb54Lrxb6gBJ1yMsExahVvvTt4fyg5Ub3wT/4Kh+2bMkuKAce6JYmKNVPttHyqZWGoXAPpdBwF1jIyzCi4uijYU7DcFqkHv7v/wo7id9D8UphpBstbyGv6idbPS//5FoeXbu6Zb4eSqEJebCQl2FERdeusNcJzfy6y9Xw0EP0LmSA4+bNyfItHTtCQ0N6QTEPpfrJ5qGsWOE8Ev+DRb5JedXiPZTu3V0JIPNQDCN8TjsNrt14GVsbutPn5ZfzP4E/5AUu7GWCUptkm2QrtdIw5B/yWrnSiU8xHkqHDtGVX1E1QTFqm7FjoUuvLiztOIBOhZSXT+3BZYJSu+QKefnDXZC/oBQ7BsUjKkH57DP3ezBBMWqVzp3dyPkP1/en4eN5+c9ZElRQwqhqHAEmKCHS1ORu1K2t7bel1vGC/AWl2DEoHlFVHE4tu2IYtcipp8KClv70mPtRcrR7UPw5FMgsKPmMhi4hFSsoInKCiNwlIn8TkcPLbQ84D0Q1fVmTMEJeYXkoURWITFd2xTBqja98BYZ9pX9hBwfNoWSbJ6OMlEVQROSPIrJERKanrD9CRGaJyGwRuTzbOVT1cVU9DzgbODVCcwOTrZ5XupBXvr28li51S28G0EKJKuRlHophOAdjr902FXZwupBXuidUE5Q2jAeO8K8QkTrgNuBIYBgwTkSGicgXROTJlJf/lnp14riyk62eV7qQV329++4E7eW1alVyRH4xWMjLMCKmubmw41JDXr16maDkQlVfAlJvu/sCs1V1jqpuBh4CjlfVd1X1mJTXEnFcDzyjqrGYuixTCfuNG13COlVQIL8S9mH1oIoq5OWJlAmKUfP84Aes3Htv9/6II9w8KUFI9VB6906WKvez557h2Bky9bl3KRkDgPm+zwuA/bLs/z3gUKCXiAxV1TvS7SQi5wPnA/Tr14/JkycXZNy6detyHjt3bldgX15+eQadOi39fP2yZZ2AA1iyZBaTJy9uc0x9/f7Mnr2SyZNn5bRh9uzh1Nd3YfLk1wr4C5KsXDmYjRt3ZNKkF+nUKY/a+Tn4z3+2B3ZhxowpLF68Je/jg7SxUTzWzqWh34gRNL75Jjz7LK889hibA8x7vc/KlWzq2JHpif/PoGXL2Km1lX898wwtXbsyBlhy0EG8t+OOEMP/YZwERdKsy3i3U9WbgZtznVRV7wTuBBg5cqSOGTOmIOMmT55MrmMXJ7Sif/898O86PZEpGj16V8aM2bXNMa6oaH/GjMmdxKuvh4EDyWlHLmbMcMsRIw76fPqFMHjlFbc86qgvtStfH4QgbWwUj7Vzafjg8cc/f3/A/vsH603TqRM9tt8++f+ZPRuA/9lzT9h+ewC2PeQQtj344LDNDYU49fJaAPhbfCCwqEy2FESmEvbpKg175Bvy6tWrcPs8oqrntWqVmwelEDExjGqjpaEh+SFoOft0ORRwPy4vbFYfJz+gLXESlFeBnUVkiIh0Ak4DniizTXnRpYvruZWaQ0lXadjDm7UxCGHmULzzhYmNkjeMJC1eN07IT1BScyjg8iieoBTbKydCytVt+EFgKrCriCwQkXNUdStwEfAsMBOYoKozQrresSJy5+qo5r31kW60fDZB6dYtv15eYdywo/RQTFAMw1GQh5I6DsXvobS0uPcx9lDKYpmqjsuw/mng6QiuNxGYOHLkyPPCPncq6QpE5gp5LViQ+7ytrW58U5iCEnbXYRMUw0hScMgrnaD4PZQYC0qcQl5VQSYPpa4uOV+On6A5lLVrXc9BC3kZRmXQJuT15S/D88/nPig1h2Ihr9omnYfijZKXNP3YggpKmIMGowx5WdkVw3C0pFYEHj8+90GZPJQKCXmZoIRMOg9l5crMN9pyCEqXLu47ayEvw4iOFq9Yn0eQkvOpOZQuXdyT6MaNFvKKC6VMymfyUNIl5MH18tqwof1A2FTCFBSR8EfL21wohtGWLb16wTPPJFfkmra3pcUlS1P3q6tz2yzkFQ9UdaKqnt8rjEEcOWhqcgLhz8FlExTvISZXaZ6w62SFXSBy40b3cGWCYhg+vvCF5PtcgrIpUVAyda4TT1As5FV7pBvcmCvkBbnDXlEISpghL6vjZRhp8I/yzRXyuvJKt0ztvVNX57wTC3nVHulK2AfxUMohKGF6KFZp2DDS4BeUXB7KnDnQ0ADnnNN2fbWHvESkW6LcvJFCasXh1lZ3sw1LUNJ1PS6EsHMoJiiGkQa/oOTyLFRh993bezLVFvISkQ4i8nUReUpElgDvA4tFZIaI/FpEdo7WzMoh1UPxKk9nCnkFnWRr1Sro0SO871LYIS8TFMNIg9+b8AQhE62t6ccWpHoolS4owD+BnYArgO1UdZD3NyJtAAAWMUlEQVSqbgv8DzANuE5EzojIxqIpdS8vSHoo2cquQNJDyVV+JeweVF7IK1fvsqCYoBhGDnLNA68KHdLckiso5BVU6g5V1XYTXKjqCuCvwF9FpGP7w+JBKUuvpIa8PC8gjJBXmDfrxkb3/dywIWlDMZigGEYOPvjAhSwy9TbN5aFUS8grVUzS5VDSCU4t0quX+054QuIJSxi9vML2ULzzhoF3nhL0zDaMyuSxx2DUqMzbg3oolS4olkMJTocO7madb8irXIISVh5l1SrXQSXIYGDDqFk+/DDztqA5lBiHvGoih1Jq/KPls1UahvKGvLzzhoHV8TKMIsnloVRAyKsmciilpqmpvYcSRi+vuIe8LH9iGGl4+GE4+eTc+9VKLy9PTETkRRHpmXj/HRH5QWJ2Rcuh+PAXiFyxwnkhmUJB3pQJ2Xp5tba6XF7cQ14mKIaRhpNO+nw++KxUQS+vfAc29lbVNSKyD3Ae0AjcFb5Z4VLKbsPQPuSVLRTUoUPuaYC9uVDCTHibh2IYJcQTg2zk8lBefNF9rnQPxccWEakHzgSuV9X/BfYI36xwKWVxSGjvoWRKyHvkKmHv6WCcQ14rV5qgGEZGggiKamZBmTMHbrgh+Tmm5Ct1NwNvA12AyxPruodqURXgeSitreEIShRjPDp2dNe1kJdhlICggpIp5OWPiVeLh6Kq9wH7AcNVdaOIDAWmRmJZBdPU5MRk7drcIS8oj6B45wvDQ7G5UAwjB7nKrkD2kJd/fYwFJZBlIiKqrkiHqq7z1qvqbOCbqfvUOv7R8kE8FG+SrUxEJShhFYhcv979XkxQDCMDWwL0WcrmofjLtsQ45BV4HIqIfE9EBvtXikgnETlERO4FzgrfvMrEXyAyriEv73xhCIqVXTGMHATp8pvNQ/Em38p1jjIT1LIjgG8BD4rIEGAVLo9SB0wCblLVt6IxsfLwPJRFi9zMjUFCXosXZ94epaAsWFD8eUxQDCMHra1umW1OlKAeShUIyndV9Sbg9sQAxj7ARlUNcUaN6BCRY4Fjhw4dWpLreR7Jf//b9nMmgnooYXdSa2yE6dOLP48JimEEJHV6Xz/ZPJQqC3ld5Hv/NVVd7ImJiPQTkSPjPFK+HN2GIVxB6d49/AeTsENeVnrFMHKQTVCyeSgVEvIKKiiDRaRH4v3vUrbdB5wKPBCaVRVOqocSRi+vKJ7+e/d2Y1w8b7xQzEMxjICE4aFUgaCsAH4pIscD9SLyZd+2/qp6NnBv2MZVKg0NbubP2bPd5zB6eUUlKKqwZk1x5zFBMYyAFJpD8T/1xTjkFVTqTgb64cqtnATcIiI3AtsBSwBU9alILKxQGhvho4/c+yAhr02bXNfbdN+VqATFX3G4mPPbXCiGEQLZPBQ/le6hqOpLqvqwqh6jqs8CpwB7Ac04kTFSaGpKdj0PEvKCzGGvKD0USD9aXhVuvBGmTMl9nlWr3N+QzZs3jJpm4kS3zDZiPpuH4nH77cmKsjEk31peAKjqLFX9kapeoKofhW1UNeCJSF0d9OyZfd9yC0q6xPyVV8Kll8Ktt+Y+j9XxMowcHHMMfPvbsHAhjB6dfp8gHsrw4dHYFxI5BUVEDhORu0Rkr8Tn86M3q/LxwlyNjem/I37KJSiZJtn6zW/guuvc93ju3NznsbIrhhEAL1Q1bVr67UE8lBjnTyBYDuUCXHmVq0WkCRfqMnLgF5RcZBOUKOZC8Ujnodx3H1xyiZvCoXt3+Pvfc5/HBMUwApArJhzEQ4lx/gSChbyWquoqVb0UOBwYFbFNVYEnJLkS8pCctTFdT69169z3rBQ5lCefhG99Cw49FO6/H3baCT75BDZuzH4eExTDCEAuMagCDyWIoHzee0tVL8eNO6koSj3BFoTnoUTZJbdnT/dAtGoVvPyym6V0773h0UfdDJPNzW6/jz/Ofh4TFMMIgHkooKp/S/l8S3TmREOpR8pDfh5KuQSlQwfX1feVV1zOcIcd4OmnoUdiCKsnKLnyKCYohhGAIB5KLkGpAg/lc7LNKW+0xROSYgXFc6qi0sLeveH5512+ZNIk6Ns3uW3IELfMJijeXChWdsUwcpDLQwkS8qp0DyWFipxTvhxUQsgLYNttna2TJsHgwW239e/vfgMfZekYHmWOxzCqinRisH59sk5XkJBXzD2UfOUudU75CSLyWgR2VTyFhLzSJeWjFpR77nH5kp12ar+tQwcXBsvmoVjZFcMISDoPpXt3GDYMZsyoCg/F5pSPiAED3Hdj0KDc+3q9vMrhoQwbln37kCEmKIYRCpnE4L334JZbqsJDsTnlI2LAAHj3XTjhhNz7durkvifZBKVcdbKam7OHvExQDCMg2byL73+/ZroNI5KUTVVdp6obE+9nq+o3U/cxHMOGBfv/i2QuYV/uOlnNzbB0afZR/GCCYhhFUwvdhhPYnPIRk01Qynmz9np6ZRqLYoJiGAHxqsVmolY8FNyc8i24OeUXich7IvIR8CEwDjen/PiIbKwJ4ioo3liUTGEvb5S9CYph5CBbpWGoCg8lkHWq+hlwOylzygPDgXGq+lZ0JtYGmSbZiougZErMlzvHYxgVg99DSSceVeChFCJ3e+C8klOBT4HdgAvDNKoWyeah9O9fens8ttvOzT6ZTVCimO/eMKqOVEFRbbu9CjyUoEn5XUTkGhF5H/gDbkrgMaq6X+J9rClHLa98iWvIS8SNRckU8iq3fYZRMfgFpaWl7Tzx4ASlwj2UoDmU94GjgZNUdR9VvV5V5ya2aebD4kE5annlS1wFBVzYK5uHYmVXDCMAuQSlpaU2PBTga8Bc4B8i8qfEE79N+Boi6QTFq5NVbkHJNrgxDvYZRkUwcGDyfTpB2bIlt4eSbnuMCDqn/GOqeiowFPg78G1ggYjcA+SY4NYIQjpBiUudrOZmWL4c1q5tv80ExTACcsEFcNBB7n1LS7KGl8fWrbk9lJiT70j59ar6gKoeA+wOTAPejcSyGiNdL6+4jPHI1tPLBMUwAlJXByee6N63thbmocScgv0nVV2hqr9X1YPDNKhW8TwUf8ePuAhKtjL2JiiGkQeeOGRKyteSh2JER7du7b9jcRGUTIMbo5zv3jCqkmyCAiYoRjikmxMlLoLSt68LyaV6KHHJ8RhGxeAXlNQcCtRuyMsIl3SC4g2bKfcNWyR91eG4CJ5hVAyeYKTLoYB5KEY4pJtkK05lTdKNRbE6XoaRJ7lCXuahGGGQbpKtuAuKeSiGkSeF5FBiPpjRjwlKTMiUQ+na1U3AVW6GDHH2eCICJiiGkTeFeCg9K2eonwlKTMgkKHG5Wacbi+IJipVeMYyAeILS2po+KZ/OQ9lmm2htChETlJhQyYISFxsNI/Z4Hkg+HkpTU7Q2hYgJSkyIu6B4gxv9Pb08Qakgj9wwykshORTzUIx8ydTLKy6C0tTk5j1J9VB69qyoTiiGUV4KyaGYh2LkS6ZeXnERFG8sSqqgxMU+w6gICsmhxKFXTkBMUGJC3ENe4MJeqSGvONlnGLGnkBxKBVHZ1gekEmZsrKuDzp2TghKXuVD8eB6KV8AybvYZRuwpJIdSQdSEoFTCjI3Qdk6U9evddy5ON+zmZjcnijdC3gTFMPLELygbN7bfnslDaW6G/v0jMyssakJQKgW/oMSxS25qTy8TFMPIE38O5bPPXPL0Jz9Jbs/kofz3v7BgQfT2FYkJSozo1i3ZyyuOgpI6FmXlynjZZxixx59D2bgRGhpgzJj229MdVwH5lfhbWEN07RpvD8UvKK2tsGZNvOwzjNiTGvJqaGhbq8tyKEZYxD3k1djoClV+9JETE1Uru2IYeZFOUPwDuUxQjLCIu6BAsqdXXO0zjFiTmkPp0qWtoFRAWCsblW19lWGCYhhVTrocinkoRhT4k/LekJm49XT2Bjfa5FqGUQC5Ql7moRhhkeqhxGUuFD/NzU70Zs92n01QDCMPLClvlIrUXl5x804g2dPrrbfc0gTFMPLAcihGqfBCXq2t8R006A1ufPNNt4yjjYYRWzzBOP98WLfOcihGdHgFIjdujK+g7LCDW779tvvu21wohpEHnnisXOmSkZZDMaLCX3E4roLSq5cbe7JhgxOTCv/+G0ZpSZ08yDwUIyr8k2zFVVAgGfaKq32GEVtSBaVLl7ZJ+Qp/Qqts66sM/yRbcRYULzEfV/sMI7akCoZ5KEZUVELIC5KCYmVXDCNPcoW8zEMxwsITlCVLYOvW+AqKhbwMo0Ash2KUCk9QFi50y7jesC3kZRgFYjkUo1SYoBhGlZMqGJ07m4diRIMnKIsWuWVcb9jNze430LdvuS0xjAoj1UOpr6+qHEp97l2MUuH18oq7h9K9O7zwAgwfXm5LDKPCyCUoFe6hmKDEiEoJeQF8+cvltsAwKpAqF5TK9q+qjIYG932qBEExDKMAUgWlrq6tiFR4yKuyra8yRFzYa8UK9zmO1YYNwyiCVMGoTwkSmYdihIkX9mpocB1ADMOoIjp2bPs5VVDMQykPIrK7iNwhIo+IyHfLbU9YeIJi4S7DqEJE4NVXk59TQ2DmoeSPiPxRRJaIyPSU9UeIyCwRmS0il2c7h6rOVNXvAKcAI6O0t5R4Pb1MUAyjSvGLiHkooTAeOMK/QkTqgNuAI4FhwDgRGSYiXxCRJ1Ne2yaOOQ54GXi+tOZHh3kohlHlZBOUCvdQytJtWFVfEpHmlNX7ArNVdQ6AiDwEHK+qvwKOyXCeJ4AnROQp4M/RWVw6PEGxhLxhVCl+QUkNeVW4hxKncSgDgPm+zwuA/TLtLCJjgK8CnYGns+x3PnA+QL9+/Zg8eXJBxq1bt67gY/Phs8+GA33YvPlTJk+eGfn14kSp2rjWsXYuDZnauWHevM9vbG+88w5rWloYk/j87vTpLK/g8EScBCWdr6eZdlbVycDkXCdV1TuBOwFGjhypY8aMKci4yZMnU+ix+TB4MEydCrvs0o8xY/pFfr04Uao2rnWsnUtDxnaePfvzt1/cd18YNerzz18YMQIq+H8TJ/9qATDI93kgsKhMtpQNy6EYRpVTxTmUOAnKq8DOIjJERDoBpwFPlNmmkmO9vAyjykmXQ/FyJxWeQylXt+EHganAriKyQETOUdWtwEXAs8BMYIKqzgjpeseKyJ2rV68O43SRYh6KYVQ56TwUb1nhHkq5enmNy7D+abIk2Iu43kRg4siRI88L+9xhY4JiGFWOP8zlve/YETZvNg/FCBcTFMOoctKFvLySLBXuoZigxAwTFMOoctKFvDxBMQ/FCJM+fdyyX231GDaM2qGKcyg1ISiVlJQ//niYMiU5b7thGFWGhbwqG1WdqKrn96qAeib19XDAAeW2wjCMyLCQl2EYhhEKmXp5gXkohmEYRh5kC3mZh2IYhmEExu+FWFLeMAzDCAXLoVQeldTLyzCMGiI15NXaWj5bQqAmBKWSenkZhlFDeB6JJyhbtpTPlhCoCUExDMOINSYohmEYRih4grJ1a3ntKBITFMMwjHLjJefNQzEMwzCKwkJehmEYRiiYoFQO1m3YMIxYc845brn//uW1o0hqQlCs27BhGLHm0ENBFYYMKbclRVETgmIYhmFEjwmKYRiGEQomKIZhGEYomKAYhmEYoWCCYhiGYYSCCYphGIYRCjUhKDYOxTAMI3pqQlBsHIphGEb01ISgGIZhGNFTX24DDMMwao7LLqv42RnTYYJiGIZRaq67rtwWRIKFvAzDMIxQMEExDMMwQsEExTAMwwgFExTDMAwjFGpCUGxgo2EYRvTUhKDYwEbDMIzoqQlBMQzDMKLHBMUwDMMIBVHVcttQMkRkKbAKSE2m9Aqwrg+wLDrrctoT1fG59i10e5A2tTYOtk8xbZxunbVzftsK+S5D5bRzrmN3UNW+gc6kqjX1Au4sZB3wWjltjOr4XPsWuj1gm1obB9inmDa2di5NO2fYpyLaudj/kf9ViyGviUWsKxXFXjuf43PtW+j2IG1qbRxsn2LaOOj1o6JW2rmcbVzs9UOzvaZCXsUgIq+p6shy21HNWBuXBmvn0lCL7VyLHkqh3FluA2oAa+PSYO1cGmqunc1DMQzDMELBPBTDMAwjFExQDMMwjFAwQTEMwzBCwQQlBERkRxH5g4g8Um5bqgkR6SYi94rIXSJyerntqVbs+xs9InJC4nv8NxE5vNz2REXNC4qI/FFElojI9JT1R4jILBGZLSKXZzuHqs5R1XOitbQ6yLO9vwo8oqrnAceV3NgKJp92tu9vYeTZxo8nvsdnA6eWwdySUPOCAowHjvCvEJE64DbgSGAYME5EhonIF0TkyZTXtqU3uaIZT8D2BgYC8xO7tZTQxmpgPMHb2SiM8eTfxlcntlcl9eU2oNyo6ksi0pyyel9gtqrOARCRh4DjVfVXwDGltbC6yKe9gQU4UXkLe/jJizzb+b3SWlcd5NPGIjITuA54RlXfKKmhJcR+pOkZQPLJGNyNbUCmnUVkGxG5A9hbRK6I2rgqJFN7Pwp8TUR+R/lLW1QDadvZvr+hkum7/D3gUOAkEflOOQwrBTXvoWRA0qzLOAJUVZcDVfslKQFp21tV1wPfLLUxVUymdrbvb3hkauObgZtLbUypMQ8lPQuAQb7PA4FFZbKlFrD2Lg3WztFT021sgpKeV4GdRWSIiHQCTgOeKLNN1Yy1d2mwdo6emm7jmhcUEXkQmArsKiILROQcVd0KXAQ8C8wEJqjqjHLaWS1Ye5cGa+fosTZujxWHNAzDMEKh5j0UwzAMIxxMUAzDMIxQMEExDMMwQsEExTAMwwgFExTDMAwjFExQDMMwjFAwQTEMwzBCwQTFMAzDCAUTFMOIGBFpEJEXE3NlICIqIn/yba8XkaUi8mQB5+4kIi+JiBV6NcqOCYphRM+3gEdV1ZskbD0wXEQaEp8PAxYWcmJV3Qw8TxXPAmhUDiYohlEkInKWiLwuIu+IyL/S7HI68LeUdc8ARyfejwMeTJyrWUTeF5F7E+d7RES6+q51ZmL92z4v5/HENQyjrJigGEYRiEgP4DJgtKruCRybsr0TsKOqzk059CHgNBHpAuwJ/Nu3bVfgzsT51gAXJM61B3AVcIiqjgAuTuw/HRgV5t9lGIVggmIYxdECNAA3ishIVV2Vsr0PkLoOVX0HaMZ5J0+nbJ6vqlMS7+8HDky8PwR4RFWXJc6xIrFsATYnxM0wyoYJimEUgapuAIYDU4A7ReSClF02Al0yHP4EcAOJcJf/tBk+S5ptHp2Bz4LYbBhRYYJiGEUgIjur6npVfQh4khTxUNWVQF0itJXKH4Gfquq7KesHi8joxPtxwMuJ988Dp4jINolrNyWW2wBLVXVLKH+UYRSICYphFMdVIjJLRN4AhgC3p9lnEsmw1eeo6gJV/W2a/WcCZ4nIO0AT8LvE/jOAXwAvisjbwG8S+x9M+7CZYZQcm2DLMCJGRPYGfqSq3wiwbzPwpKoOz+P8jwJXqOqsgo00jBAwD8UwIkZV3wT+6Q1sDJNEL7LHTUyMOGAeimEYhhEK5qEYhmEYoWCCYhiGYYSCCYphGIYRCiYohmEYRiiYoBiGYRihYIJiGIZhhIIJimEYhhEKJiiGYRhGKPx/f5+ROLW9zDUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Max diff for s = 0.1-100: 0.04527382016672736\n",
      "Max diff for s = 50-250: 0.038796943897711024\n",
      "BAO peak height = max (s^2 xi) for s = 140-160 Mpc: 51.71152471864682\n",
      "Average Delta(r^2 xi) for s = 140-150 Mpc: 0.010881527323687957\n"
     ]
    }
   ],
   "source": [
    "frac_diff1 = []    \n",
    "frac_diff2 = []\n",
    "abs_diff1 = []    \n",
    "abs_diff2 = []\n",
    "\n",
    "s1 = rsd_bench[0:40,0]\n",
    "s2 = rsd_bench[40:140,0]\n",
    "\n",
    "xi_ccl1 = ccl.correlation_3dRsd_avgmu(cosmo,a,s1,beta)\n",
    "xi_cosmomad1 = rsd_bench[0:40,1]\n",
    "\n",
    "xi_ccl2 = ccl.correlation_3dRsd_avgmu(cosmo,a,s2,beta)\n",
    "xi_cosmomad2 = rsd_bench[40:140,1]\n",
    "\n",
    "for i in range(len(s1)):\n",
    "    frac_diff1.append(np.abs(xi_ccl1[i]/xi_cosmomad1[i] - 1.))\n",
    "    abs_diff1.append(np.abs(s1[i]*s1[i]*(xi_ccl1[i]-xi_cosmomad1[i])))\n",
    "    \n",
    "for i in range(len(s2)):\n",
    "    frac_diff2.append(np.abs(xi_ccl2[i]/xi_cosmomad2[i] - 1.))\n",
    "    abs_diff2.append(np.abs(s2[i]*s2[i]*(xi_ccl2[i]-xi_cosmomad2[i])))\n",
    "    \n",
    "r2xi = np.array(s2*s2*xi_cosmomad2)\n",
    "\n",
    "# plot relative and absolute difference in r^2 xi(r)\n",
    "plt.plot(s1, frac_diff1, 'b-')\n",
    "plt.plot(s2, frac_diff2, 'r-')\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.xlabel(r'$s$ (Mpc)')\n",
    "plt.ylabel(r'$\\Delta \\xi(s) / \\xi(s)$')\n",
    "plt.grid(which='minor')\n",
    "plt.title('Relative difference') \n",
    "plt.grid(which='both')\n",
    "plt.savefig('benchmark_abs.pdf',bbox_inches = 'tight')\n",
    "plt.show()\n",
    "\n",
    "plt.plot(s1, abs_diff1, 'b-')\n",
    "plt.plot(s2, abs_diff2, 'r-')\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.xlabel(r'$s$ (Mpc)')\n",
    "plt.ylabel(r'$\\Delta (s^2 \\xi(s)) $')\n",
    "plt.grid(which='minor')\n",
    "plt.title('Absolute difference') \n",
    "plt.grid(which='both')\n",
    "plt.savefig('benchmark_abs.pdf',bbox_inches = 'tight')\n",
    "plt.show()\n",
    "\n",
    "#print abs_diff\n",
    "print(f'Max diff for s = 0.1-100: {np.amax(abs_diff1)}')\n",
    "print(f'Max diff for s = 50-250: {np.amax(abs_diff2)}')\n",
    "apex = np.amax(r2xi[(s2<160)&(s2>140)])\n",
    "print(f'BAO peak height = max (s^2 xi) for s = 140-160 Mpc: {apex}')\n",
    "# find and print the average of Delta(r^2 xi) in the BAO peak region\n",
    "avg_value = np.average(np.array(abs_diff2)[(140<s2) & (s2<150)])\n",
    "print(f'Average Delta(r^2 xi) for s = 140-150 Mpc: {avg_value}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
